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Quark mass thresholds in QCD thermodynamics 
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We discuss radiative corrections to how quark mass thresholds are crossed, as a function of the 
temperature, in basic thermodynamic observables such as the pressure, the energy and entropy 
densities, and the heat capacity of high temperature QCD. The indication from leading order that 
the charm quark plays a visible role at surprisingly low temperatures, is confirmed. We also sketch 
a way to obtain phenomenological estimates relevant for generic expansion rate computations at 
temperatures between the QCD and electroweak scales, pointing out where improvements over the 
current knowledge are particularly welcome. 



PACS numbers: ll.10.Wx, 11.15.Bt, 12.38.Bx, 98.80.Cq 



o ■ 
o : 

(N . 



I. INTRODUCTION 



Besides being of fundamental theoretical interest to finite temperature field theory, the thermodynamic pressure of 
the Standard Model, as a function of the temperature T and of various chemical potentials /i,, has several potential 

. phenomenological applications. Most notably it dictates, through the Einstein equations, the expansion rate of the 
radiation dominated Early Universe. The expansion rate in turn determines when various dark matter candidates 
decouple, thus fixing their relic densities: fine details of the pressure could become observable for instance if dark 
matter is made of electroweak scale WIMPs 00 or of keV-scale sterile neutrinos 00. Furthermore, the pressure is 
in principle visible in the present-day spectrum of the gravitational wave background that was generated during the 

i inflationary epoch [5j . More generally, the pressure incorporates the fact that the Standard Model possesses a trace 
£C) ' anomaly, i.e. ^ 0, which in turn can influence many kinds of gravity-related cosmological scenarios (for recent 

, examples, see Refs. 0). 

Apart from cosmology, the pressure is potentially also relevant for the hydrodynamic expansion that the dense mat- 
ter generated in current and upcoming heavy ion collision experiments may undergo. In this case there is some room 
for caution, however, since the issue of whether local thermodynamic equilibrium is reached remains controversial 0. 

Qh| Given that the biggest theoretical challenges are related to strongly interacting particles, considerable efforts have 
been devoted to the determination of the QCD part of the pressure over a course of years. Denoting by g the renor- 

qj ' malised strong coupling constant, perturbative corrections to the non- interacting Stefan-Boltzmann law have been 
^ ! determined at relative orders 0(g 2 ) 0, 0(g 3 ) 0, 0(g 4 ln(l/ 5 )) 0, 0(g 4 ) [ll|, 0(g 5 ) [H, and 0(g 6 ln(l/ 5 )) 0, as 
a function of the number of colours, N c , and the number of massless quark flavours, Nf. The first presently unknown 
order, 0(g 6 ), contains non-perturbative coefficients 0,0, but those can also be attacked 0,^). All orders of g 
are available in the formal limit of large Nf • These results have been extended to the case of finite quark chemical 
potentials [Hl20ll2l|. and a similar computation has recently also been finalised for the weakly interacting part of 
the Standard Model, at temperatures higher than the electroweak scal e l22|- Moreover, the fact that several orders 
are available allows to experiment with various kinds of resummations |23l |24| . 

Surprisingly, however, relatively little seems to be known about the dependence of the QCD pressure on the quark 
masses m*, i = l,...,Nf. While the non-interacting Stefan-Boltzmann law is readily extended to this situation, it 
in fact appears that even the first non-trivial term, 0(g 2 ), has not been exhaustively investigated in the literature 
(see, however, Ref. 25] ). In principle this term has of course been available since almost 30 years |9|, but in explicit 
form only in a renormalization scheme for quark masses which differs from the current standard, the MS scheme. 
Furthermore, no general numerical evaluation of the basic integrals appearing has been presented, as far as we know. 
For T = but fa ^ 0, the full 0(g 2 ) analysis has also only been carried out recently |26j . 

Several probable reasons for the apparent lack of interest can surely be envisaged. First of all, the dependence on 
Nf is known to a high order in the massless case, and interpolating between integer values of Nf should give much of 
the information that we may need for the massive case. Second, including quark masses turns out to be technically 
cumbersome 



Third, there are several indications, for instance from considerations of the baryon chemical poten- 
tial 0, [2(j and of mesonic correlation lengths [2^ , that the convergence is much better in the quark sector than in 
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the gluonic sector, so that the lowest non-trivial order may already provide sufficient accuracy. Nevertheless, we feel 
that the last assumption deserves to be checked, at least at the next-to-leading order (NLO) 0(g 2 ), and this is the 
purpose of the present paper. 

In short, our general philosophy will then be to account for the gluonic contributions to the highest order available, 
0(g 6 ln(l/<7)), and consider the change that quarks with finite physical masses inflict on this result at NLO, 0(g 2 ). 
We do find that the quark mass effects at NLO are not too different from those at the leading order, O(g ), such 
that the philosophy of terminating at 0(g 2 ) is at least self-consistent. Nevertheless, we also outline the procedure for 
determining the quark mass dependence up to the order 0(g 6 ). 

Apart from the theoretical goals mentioned, we also wish to sketch certain phenomenological results in this paper. 
Consider the temperature evolution of an expanding system in the case of cosmology, for instance. Einstein equations 
then lead to (see, e.g., Ref. |2q) 

ldT_ l 24ne(T) s(T) 

T dt V m%y c(T) ' { ' 

where t is the time; we assumed the Universe to be flat (k = 0); and we ignored the cosmological constant. All the 
quantities appearing here follow from the pressure: s(T) — p'(T) is the entropy density, e(T) = Ts(T) — p(T) is the 
energy density, and c(T) = e'(T) = Tp"(T) is the heat capacity. We wish to present our favoured "fits" for all these 
functions for temperatures between the QCD and electroweak scales, indicating where further work is required. 

The plan of this paper is the following. We start by elaborating on the basic formalism in Sec. [HI discuss quark 
mass thresholds in Sec. IIIII present a phenomenological evaluation of the various thermodynamic functions relevant 
for physical QCD in Sec. IIVI include weakly interacting particles in Sec. [VJ an d conclude in Sec. IVII 



II. BASIC FORMALISM 



In order to determine the basic thermodynamic quantities of the Standard Model, all of which can be derived from 
minus the grand canonical free energy density, or the pressure p(T, fj,), where fj, collects together the various chemical 
potentials associated with conserved global charges, 1 we make use of the framework of dimensionally reduced effective 
field theories This framework allows to organise the computation in a transparent way, and implements 

various resummations of higher order effects. We start by briefly reviewing certain aspects of the general framework 
for QCD; further details can be found in Ref. [l3j . 

Dimensional reduction proceeds by first integrating out the "hard modes" , with momenta or Matsubara frequencies 
of order 2nT . This produces an effective theory |29j |. called EQCD |3l]], which is a three-dimensional SU(A r c ) gauge 
theory with a scalar field in the adjoint representation. The effective theory has a certain number of couplings, 
parametrised by f unctions denoted by a^\...a^ [contributing up to 0(g e ln(l/g))] and /3ei---/3e5 [contributing at 
0(g 6 )] in Ref. [13j: in the following we explicitly specify the definitions for only a subset of them. These parameters 
contain all the information concerning the hard modes. Assuming the use of dimensional regularization, we denote 
by a M f ■ /3eT parameters from which the f / e-divergences have been removed by the MS-prescription. 

To proceed, we need to specify explicitly the effective mass parameter m\ and the effective gauge coupling g\ of 
EQCD at NLO: 

2 4 

- 2 _ ^3 _ 2 MS , 9 MS (t> \ 
3 — rp2 — 9 a Ei + (4^)2 a E6 ) \ z ) 

2 4 

f-2 _ 93 _ 2 , _9 ms / o\ 

.93 - T - 9 + (47r )2 a E7 ■ W 

Both parameters are renormalization group invariant up to the order computed, i.e., the dependence on the scale 
parameter p, is of order 0(g 6 ). 

With this notation, the physical pressure of hot QCD can be written in the form 

PQCD = Phard + Psoft , (4) 



1 The notation p(T, (jl) always implicitly refers to the ultraviolet finite difference p(T, /x) — p(0, 0). 
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where phard represents the contribution of the hard modes (by definition containing both all direct hard contributions 
to the pressure, and all finite terms emerging from products like e • 1/e), while p so h represents the contribution of the 
soft modes. Up to the accuracy 0(g 6 ), phard can conveniently be expressed as 



Phard 



MS i ~2 MS 
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(Air) 4 
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4ttT 



*hard 



(5) 



where d A = N^ — 1, C A = N c , and we have separated a term on the last line which cancels the /x-dependence of p so ft 
at 0(g 6 ). The function A hard , 
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The function A so f t reads 
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where /3m can be found in Ref. 32], and a numerical estimate of (3q in Ref. |l7j . 

Let us stress that the formulae presented apply independently of whether quark masses are included or not: all 
quark mass effects can be incorporated in the perturbative functions a E f...a E7 , /3 E4 .../3 E f . In particular, the non- 
perturbative numerical value 0q and the contribution from the Debye scale /3m in Eq. JSJ are "universal" . 

In the following, we refer to the various orders of the weak-coupling expansion according to the power of <?3,?7i3 
that appear, with the rule 0(m 3 ) = 0(53) = 0(g)- In other words, "0(<?")" denotes 0(<? 3 -fc m 3 ) in the expression 
constituted by Eqs. 10} , ©, 10. If g%,rh\ were to be re-expanded in terms of g 2 , the result of Ref. pjj would be 
reproduced up to 0(g 6 ). In practice, however, it is advisable to keep the result in an unexpanded form, because this 
makes it more manageable, and because the unexpanded form introduces resummations of higher order contributions. 



III. QUARK MASS THRESHOLDS IN THE PRESSURE 

In the absence of an explicit 0(g 6 ln(l/g))-computation with ^ 0, we are forced to estimate the effects of finite 
quark masses by other means. In order to do this, we take as the starting point the limit Nf = 0, where all quarks are 
treated as infinitely heavy. Subsequently, the masses of a number Nf of them are lowered to their physical values. In 
this process the pressure at any given temperature increases. The goal is then to estimate the "correction factors" that 
describe this increase. Note that this philosophy corresponds to the "unquenching" of quark effects that is regularly 
attempted in lattice simulations. 

The philosophy just suggested is of course not unique. For instance, one could also start from the "opposite" limit, 
taking as the starting point the pressure for a fixed Nf, but with vanishing quark masses, and then increase the masses. 
In fact, one could perhaps even define an effective non-integer massless Nf by evaluating the massive a E f (Eq. (jl L|> 
below), fitting it to the massless formula (Eq. (A.l) in Ref. [L|)j and using the resulting Nf in the massless result of 
0(g 6 ln(l/<7)) |33(. We prefer, however, to take here the infinitely massive case Nf = as the reference point, since 
this limit can be treated with more confidence than the chiral limit, as we describe in the next section, and since this 
limit offers us a convenient way to probe the convergence of our recipe, as we now explain. 

The simplest thinkable way to estimate the correction factors due to non-infinite quark masses would be to multiply 
the Nf = result with the change indicated by the Stefan-Boltzmann law, i.e. by a E f (JVf)/a E f (0) 0. The first non- 
trivial improvement of this philosophy is then to determine the functions a E ^ , a E 2 1 c^e? m the general massive case, 
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FIG. 1: Left: the pressure for Nf = 0, 3, 4, at O(g ) and (D(g 2 ). Right: the "correction factors" accounting for the effects 
of quarks, at O(g ) and 0(g 2 ) (cf. Eq. 1101 1. They grey bands indicate the effect of MS scheme scale variations by a factor 
0.5 ... 2.0 around the "optimal" value. It is observed that while the 0(g 2 ) corrections are of order 20. ..30% in the pressure, 
they are of order 5% in the correction factors for 7Vf = 3, and even less for the physical case Ni = 4. 



which allows us to evaluate the order 0{g 2 ) result for the pressure, viz. 
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Afterwards, we may modify the Nf = result with a correction factor, 
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Comparing the outcome of this 0(g 2 ) recipe with the corresponding 0(g ) recipe allows to probe the convergence. 
Note that it is important to also determine ajf| , since only this way can the renormalisation scale that appears in g 2 
be reasonably fixed (cf. Eq. J3J). 

We thus proceed to compute a^f, a^. We do this in full generality, keeping Nf, N c , the quark masses rrii, 
and the chemical potentials fii as free parameters. The quark masses and the strong gauge coupling are renormalised 
in the MS scheme. Some details concerning the computation are collected in Appendix A. As final results, we obtain 
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(12) 
(13) 



where the functions F\, ■■■,F i and some of their properties are detailed in Appendix B. 

To estimate the numerical importance of the 0(g 2 ) corrections, we need to assign a value to all the parameters that 
appear. Following a simple-minded logic, we use 1-loop running, 

9C F 

247T 2 ,_ Jln^ef/Ajjs)" 



{HCa 

n _ 



47» ln(/2/A^) ' 



rriiin) = mi(^ r ef) 



M/V A ms) 



where Tp = Nf/2, CV = (iV,r — 1)/2A C , p, Te f = 2 GeV. The quark masses at p, = /i ro f are taken from Ref. 
To choose /z, we apply the principle of minimal sensitivity criterion for the parameter g 2 , as suggested in Ref. 3 
Furthermore, for illustration, we set A^- = 200 MeV. 
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The outcome of this procedure is shown in Fig. ^ f° r Mi — 0- It is observed that while the 0(g 2 ) corrections are 
of order 20. ..30% in the pressure (left panel), the "correction factors", i.e. the ratios in Eq. (JTUJ, only contain 0(g 2 ) 
corrections of order 5% for Nf = 3, and even less for the physical case Nf = 4 (right panel). This implies that the 
quark mass dependence of the pressure probably converges faster than the weak-coupling expansion as a whole. 

Finally, we note from Fig. fright) that the charm quark contribution starts to be visible already at fairly low 
temperatures. At leading order, the quark mass dependence is determined by the function F\ (cf. Eq. (|11J) ). which at 
low temperatures has the familiar classical form 

K&oH&y -(-?)■ « i5 > 

It is observed that F\ obtains 5% of its asymptotic value 77r 2 /720 at temperatures as low as T k m/5. (For the 
precise numerical values of F\, see Fig. |3J) As Fig. fright) shows, the onset of a visible charm mass dependence is 
postponed to about T ~ 350 MeV at 0(g 2 ), but the basic pattern remains unchanged. 

IV. PHENOMENOLOGICAL RESULTS FOR QCD 

We now move from fairly well-defined expressions towards phenomenology. The goal is to present, where possible, 
an educated numerical guess for the physical QCD pressure. We set all chemical potentials to zero in the following. 

The general philosophy we adopt is that, for temperatures above the deconfinement transition, the weak-coupling 
expansion needs to be evaluated up to the order where the dominant contributions from all the different scales 
(2ttT, gT, Q 2 T) have made their entrances. There is some support for such a recipe from a number of non-trivial 
observables [3(||33. ln practice, this means that the QCD pressure would need to be evaluated up to 0(g 6 ). 

Unfortunately, one of the 0(g 6 ) terms, parametrised by (3^f in Sec.|nJ remains completely unknown for the moment 
even in the massless limit: it is a function of Nf and requires a perturbative 4-loop computation. This introduces a 
certain unknown "constant" into the prediction. We propose to fix the constant by the following recipe. 

Let us start by considering the case Nf = 0, N c = 3. Then the expressions in Sec. ^depend on only two parameters: 
on p,/T (through a^f...a^, ^ei— /?es)' an( i 011 mMms (through g 2 (p,)). It so happens that the dependence on p,, 
which formally cancels up to the order of the computation, is numerically non- monotonous (see, e.g., Ref. [2j|), so 
that the specific choice is not terribly important, as long as we are close to the extremum. In practice we choose JijT 
according to the principle of minimal sensitivity criterion for the parameter g 2 , as already mentioned. Thereby the 
results only depend on T/A^g- and on the unknown 0(g 6 ) terms, contained in Ahard and A so ft, defined in Eqs. ©, 10. 
It is important to note that once jl/T has been fixed, the A's can be treated as temperature-independent constants. 
It is furthermore convenient to combine them into a single term, 2 

A hard + d A C A A soit ee d A C A A . (16) 

In order to now eliminate the dependence on A, we "match" the perturbative prediction to 4d lattice simulation 
results for the case Nf = 0, where the continuum limit has been reached with reasonable precision [3^.l3^|. It should 
be stressed that this step is purely phenomenological: in principle A is computable from the theory. On the other 
hand, there is every reason to expect that results obtained through the dimensionally reduced framework do match 
4d lattice results as soon as T>2T C , where T c is the temperature of the deconfinement phase transition (see, e.g., 
Refs. [3(|, Moreover, that a family of functions specified by a single parameter should match a given 

function for a whole range of argument values, provides for a non-trivial consistency check. 

Now, lattice results are usually presented in terms of T/T c , rather than T/A^g. We thus need a value for T c /A^; 
we use T c /A^3s ss 1.20 which appears to be consistent with all independent determinations (cf. Ref. Sec. 4.2). 
After this choice, an excellent match can be obtained (we do this by minimising the difference squared of the function 
values in the range T > 3T C ), with a value A w —3.287 (cf. Fig. 0). In the following we will take the cubic spline 
interpolation shown in Fig. [3 as the "starting point", which will then be "corrected" by the effects of quarks. 

To now include quarks, we simply multiply the result just obtained by the correction factor in Eq. 110(1 . We should 
expect this construction to work the better the higher the temperature, but surely at least T > 200 MeV is required. 

It needs to be noted, however, that like in Fig. ^ the evaluation of the correction factor necessitates fixing A^g 
in physical units. This exercise is non-trivial. We again choose a purely phenomenological but rather convenient 



2 We note that A differs by a certain constant from a similar constant employed in the figures of Ref. Il3l . 
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FIG. 2: A phenomenological interpolating curve (solid line) for the QCD pressure at Nf = 0. In the perturbative curve (grey 
band) the unknown 0(g 6 ) constant has been adjusted so that lattice data (closed squares [38|l is matched, once T > 3.6A^Yg. 



procedure, which makes use of the pressure produced by the full set of hadronic resonances |34j . Indeed, it has been 
demonstrated recently that if the resonance masses are tuned to correspond to the quark masses accessible to current 
lattice simulations, the resulting "resonance gas pressure" works surprisingly well even for temperatures deep into the 
crossover region [43j . 

Thus, we tune such that our analytic recipe and the first derivative thereof match the resonance gas result (in 
which the temperature is automatically measured in physical units) at a certain temperature. Examples are shown 
in Fig. Ofleft). The values of A^ that result depend slightly on Nf and on variations of quark masses within their 
experimental errors, but the typical range is A„g w 175. ..180 MeV. We should stress that this matching is of course 
rather arbitrary, but it does produce qualitatively reasonable results with, for instance, an inflection point on the side 
T > T c ~ 175 MeV, as suggested by lattice results |4^-|47|. 

Naturally, the resonance gas results cannot really be trusted in quantitative detail for temperatures above, say, 
150 MeV. Therefore, for a certain interval (which we choose to be T — 150. ..350 MeV, and shade in all figures), the 
results remain to be established by lattice simulations. The matter becomes even more urgent, when one considers 
derivatives of the pressure, to which we now turn. 

Apart from the pressure, its first and second derivatives play an important role, as already mentioned in connection 
with Eq. HJ. There are various ways of presenting the information contained in these derivatives: we may for instance 
parametrise the physical observables e(T), s(T), c(T) through effective numbers of bosonic degrees of freedom, 
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or we can consider dimensionless ratios like 
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Both the "equation-of-state" w(T) and the sound speed squared c^(T) equal 1/3 in the non-interacting limit. The 
deviation of the parameter w(T) from 1/3 is proportional to the trace anomaly, sometimes also called the interaction 
measure. 

Results for all of these quantities, based on our interpolation, are shown in Fig. EH middle, right). It can be seen 
that quantities involving derivatives show a significant amount of structure around the QCD crossover, even if there 
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FIG. 3: Left: Phenomenological interpolating curves (solid lines) for the QCD pressure at iVf = 3,4. The shaded interval 
corresponds to the transition region where the results can be reliably determined with lattice simulations only. Middle: 
Seffi ^effi *eff as defined in Eq. $17) . for Nf — 4. Right: the equation-of-state parameter w and the speed of sound squared (? s , 
for the same system. 



were no singularities. We remark that to smooth the behaviour we have evaluated p(T, fi) with a relatively sparse 
temperature grid in the critical region. 

Clearly, it is important to correct the results in the "shaded re gion " by using results from future lattice simulations 
of the type in Refs. I n particular, the recent Refs. [irj l47j display direct results for c 2 and w, respectively. 

Unfortunately, it does not appear that these results would be useful for our present purposes: they for instance fail 
to reproduce the significant rise in w and c 2 that is seen in Fig. [Sfright) at temperatures down from the critical one, 
displaying rather a much deeper dip (down to ~ 0.1) around the critical region, and then rising at most slightly as the 
temperature is lowered. Therefore, it could be feared that the dip itself is affected by the unphysically heavy quark 
masses that are used in the simulations. 

We finally comment on the peak visible in i e g in Fig. [21 middle). While the details are of course not captured by 
our phenomenological recipe, the fact that a peak exists in the heat capacity is not unexpected for rapid crossovers; 
in second order phase transition, the heat capacity even diverges as T — > T c . 



V. PHENOMENOLOGICAL RESULTS FOR THE STANDARD MODEL 



While in heavy ion collisions at most strongly interacting particles have time to thermalise, the expansion rate is 
much smaller in cosmology, so that all Standard Model degrees of freedom do reach thermal equilibrium, and remain 
thermalised until neutrino decoupling at around T ~ MeV. Therefore, their contributions need to be added to the 
QCD pressure. In practice, we count gluons and the four lightest quarks as the QCD degrees of freedom, while the 
bottom and top quark are treated as part of the "weakly interacting" sector, so that the result splits into a sum of 
two terms. 

We will assume that it is sufficient to treat the weakly interacting sector at 1-loop level. That is, we construct 
the free energy density / in the presence of a Higgs expectation value v, temperature T, and chemical potentials /ij, 
according to 

/(«,!» = ~v 2 (jl)v 2 + ^\(fi)v i + J2<riJi(m(v),T,»i) , (21) 

i 

where the sum extends over all physical degrees of freedom, with their proper degeneracies; <7i = +1 (—1) for bosons 
(fermions); and the tree-level masses rrii(v) depend on v in the standard way (it is sufficient at this order to work in 
unitary gauge). For scalar (J s ), vectors {J v ) and fermions (Jf), 

Js = £te**l*(l-e-^)^, (22) 
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FIG. 4: Left: The Standard Model pressure for ran = 150 GeV, 200 GeV. The shaded intervals correspond to the QCD and 
electroweak transition regions. Middle: g e g, h e g, i e g as defined in Eq. 1171 . for van = 150 GeV. Right: the equation-of-state 
parameter w and the speed of sound squared c? s , for the same system. Various sources of uncertainties are discussed in the text. 



The renormalised pressure is then given by 

p(T, (j,) = min t J(w, 0, 0) - min t ,/(w, T, (i) . (25) 

The renormalised pressure depends on a number of parameters denned in the MS scheme: the Higgs potential 
parameters ^ 2 (/i), A(/2) (cf. Eq. and the weak gauge and the top and bottom Yukawa couplings g^ift), h^(p,), 

hl(fl) (through the tree- level masses). The first four of these we express through the Fermi constant and the W , 
Higgs, and top pole masses, employing the explicit 1-loop relations listed in Ref. [3(J, while the last one is fixed 
through the bottom mass in the MS scheme 34] . Given that the electroweak theory contains a multitude of scales, 
both zero temperature and thermal, we simply choose a fixed ft = 100 GeV for the weakly interacting part of the 
pressure (we have varied the scale by a factor 0.5 ... 2.0, and seen that the dependence is invisible on our resolution). 

Let us remark that Eq. I|25|) suffers from the problem that it leads to a first order electroweak phase transition 
at a certain temperature, while there is none if the theory is treated more carefully psl |49|. In practice this does 
not lead to any serious complications, however: we again smooth the behaviour by evaluating p{T, /x) with a sparse 
temperature grid around the critical region. In our figures, we shade the corresponding temperature interval, where 
our estimates are qualitative at best. 

With these reservations, the whole Standard Model pressure, and the parameters defined in Eqs. l|17fl . (|19|) . I|2U|I . 
are shown in Fig. 

At temperatures above the electroweak scale, our results are already v ery close to the ideal gas results. Recently, 
higher loop corrections in this region have been considered in some detail |22| . The authors find a rather more signifi- 
cant deviation from the ideal gas value, due for instance to the top Yukawa coupling. We have not implemented these 
corrections, however, since they would require a correspondingly higher order computation in the broken symmetry 
phase. Though such a computation exists in principle up to 2-loop level j50j, we did not consider the non-trivial 
challenges posed by its numerical evaluation for general masses to be worth tackling at present, given that in the 
quantities plotted in Fig. fright), the 2-loop contributions (which do not contribute to the trace anomaly on the 
symmetric phase side) are expected to largely cancel out. Nevertheless, it would be important to finalise this compu- 
tation, if physics is made with the temperature interval T = 10... 100 GeV, where the W , Z° bosons and top quark 
cross their mass thresholds. 

Finally, once a definite Higgs model is available, it will of course be important to carry out lattice simulations for 
the transition r egio n. Fortunately, for the electroweak theory this can be achieved within the dimensionally reduced 
effective theory |30j . whereby also all fcrmions with their physical Yukawa couplings can be fully accounted for. 



VI. CONCLUSIONS 



The functional dependence of the QCD pressure at a high temperature T on most of the parameters of the theory 
(number of colours N c , number of active massless flavours Nf, and quark chemical potentials jii) is known up to relative 
order (D(g 6 lri(l/g)), while the dependence on the quark masses rrii has gained much less interest. The purpose of 
this note has been to study whether it indeed is justified to consider the effects of finite non-zero quark masses at the 
level of the non-interacting Stcfan-Boltzmann law (i.e. 0(g )), as has been the standard procedure. For this purpose, 
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we have determined the corrections of order 0(g 2 ) in full generality in the MS scheme, and presented a numerical 
evaluation of all the integrals that appear in this result. 

We find that while the 0(g 2 ) corrections are in general 20. ..30% (Fig.EJleft)), they are numerically at most 5% for 
the quark mass dependence (Fig. fright)). This is perhaps in accord with previous observations according to which 
quarks are fairly perturbative as soon as they are deconfined, even though gluons do display strong interactions up 
to very high temperatures. 

Finally, we have sketched educated "guesses" for the thermodynamic quantities that play a role in various physical 
contexts, for temperatures between the QCD and electroweak scales. For the case of heavy ion collisions, in particular, 
it is perhaps relevant to keep in mind that if the charm quark does thermalise, it has a rather significant effect even 
at relatively low temperatures (Fig. EJleft)). Of course, it is by no means clear whether such a thermalization should 
take place in practice |5l| . 

In order to improve on our QCD results, the missing perturbative 0(g 6 ) computations and, naturally, lattice 
simulations in the transition region, with physical values of the quark masses, remain to be completed. 

For the full Standard Model, we have presented similar guesses for the various quantities that are relevant for expan- 
sion rate and particle decoupling computations (Figs. 0J. Although the deviations from previous phenomenological 
estimates that have appeared in the literature 0, Q are in general fairly small, we nevertheless hope that our results 
help for their part to gauge the systematic uncertainties that still exist in these quantities. 

In particular, we have stressed the need for repeating the computations of Ref. [22j in the broken symmetry phase 
and, of course, the need for effective theory lattice simulations in the transition region, once the electroweak model / 
Higgs mass is known. 
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APPENDIX A: OUTLINE OF THE COMPUTATION 

In this Appendix we present a few details for the computation leading to Eqs. (|ll[l - i|13[) . We concentrate on the 
fermionic contributions; the bosonic ones are elementary. 

We write the fermionic contributions in terms of the renormalised gauge coupling g 2 and the renormalised quark 
masses rrij , i = 1 , . .. , Nf . The master integrals emerging are 

1 



h = Xp2> ( A1 ) 



J f (m, M ) = It ln(P f 2 + m 2 ) , (A2) 



2 



7 f (m, M ) ee 4 1 (A3) 
7p f P f + m 2 

HAm, u) = T — ^ = = — , (A4) 

T P{ ,Q f {Pf + m 2 )(Q 2 + m 2 )(P f - Q f ) 2 

where P^, Pf denote bosonic and fermionic Matsubara four-momenta, respectively, and Pf = Pf + (— iji, 0) includes 
the chemical potential /i. With this notation, the fermionic contributions to the parameters of interest are 

Nt 

T 4 4i = ACa^Mtth,*), (A5) 
»=i 

N f 

T 4 a f E2 = 2C A 5 1 m 2 ^2m 2 I i (m i ,fi i ) + 

i=l 

Nt ( d — 1 1 
+ cUyN 2 2I b - I t {muHi) If(m l ,fii) + 2m 2 H f (m i ,n i )'>, (A6) 

»=i ^ ' 
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(47T) 



2 "E7 



= <W 



dl{(mi,fj,i) 



3 ^ 

i=l 



(A7) 



where Sim 2 , 8\g 2 are counterterms defined by writing the bare mass parameter and gauge coupling as 



mf(l- 



g 2 8\m 2 ), g 2 B — g 2 (l + g 2 S±g 2 ); it is understood that only the fermionic part of 5±g 2 is considered; and d = 3 — 2e. 

The next step is the evaluation of the Matsubara sums appearing in the master integrals. For the 1-loop structures 
(Eqs. (|A1I) (|A3|) ). it is straightforward to obtain 



h = 
J f (m, p) = 
If(m,/x) = 



d d p 1 



{2ttY 


2|P| 


d d p 


-P 2 


(27r) d 


2dE 


d d p 


1 r 



l + 2n B (|p|)J , 
1 — n F (E — /i) — n F (E 
1 - n F (E - fi) - n F (E + 



(2n) d 2E 

where we have carried out a partial integration after the sum in Jf, and 

1 ,„ 1 



M) 



n B (E) 



n F (E) 



eP E + l 



(A8) 
(A9) 
(A10) 

(All) 



As is well known, the momentum integral in Eq. (|A8() can be carried out explicitly, I b = n d / 2 T d r(l-d/2)({2-d)/2?r 2 T, 
but the ones in Eqs. JX9), (|AT0)l with m, /i ^ cannot in general be integrated in closed form. 

For the genuine 2-loop integral Hf in Eq. I|A4|I the sums are slightly more com plic ated, so we give here some 



details. The method we employ follows the standard procedure |9( (see also Refs. [52t |53| 1. The twofold sum over the 
Matsubara modes is first written as a threefold sum with a Kronecker delta-function, and the delta-function is then 
written as S(po) = T J^dx exp(ipox). The sums can now be performed: 



T 



e ip h x 



E 2 

e ip f x 



^P 2 

Pf 1 



E 2 



1 

2E 
1 

2E 



n B (E) 



rE 



n F (E + n)e^-^ E+l3ti - n F (E - /i)e 



.rE 



(A12) 
(A13) 



where pb, Pf denote the bosonic and fermionic Matsubara frequencies, respectively; and fit = pi — ifi. The integral over 
x is then simple. All the exponents appearing can be written in terms of inverses of the distribution functions tib, 
n F , and multiplying with their explicit appearances, we are left with at most quadratic products of the distribution 
functions, and fractions containing the three "energies" . 

The fractions containing the energies can be organized in a transparent form, once we introduce the zero-temperature 
objects 



it I 2 2 2\ 

U(Q 2 - : mlm 2 2 ) - 
A(Q 2 -m 2 ) = 



'P 



(27T) 4 - 2e 
d 4 ~ 2£ P 



j4-2£ 



1 



(2tt) 4 - 2£ [P 2 
1 



m 2 ][Q 2 + m 2 ][{P - Q) 2 + m 2 ] 



(27r) 4 - 2£ [P 2 + m 2 ][{P - Q) 2 + ml 
1 



Q 2 



(A14) 
(A15) 
(A16) 



Indeed, carrying out the integrals over Po,Qq in these functions, one obtains similar energy fractions. Making 
furthermore use of the 0(4 — 2e) rotational invariance of the Q-dependence in Eq. (|A15|I . which is present once also 
the integration over p is performed, the various fractions can be identified with each other. 

In order to write the subsequent result in a compact but generic form, we introduce the notation 



pf + m 2 , P i = {E i ,p i ) 



Pi ■ Pj 



and denote 



n±(Ei) = 



_ / n B {Ei) 



for bosons (= E 3 ) 



-n F (Ei ± n) for fermions (= E\, E 2 ) 



(A17) 



(A18) 
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Then, allowing for generality for three different masses, like in Eq. I|A14(1 . 



d d Pl n a {Ei) 



E E 
E E 



n(-m^;m 3 2 ,m 2 fc ) + 



(2ir) d 2Ei 
d d Pl fd d p 3 ■ n a (Ei)n T (Ej) 



:2n) 



^EiEj 



^[-(aPi-rPj^ml], 



(A19) 



where X^j^fc = S(j j fe)=(i 2 3) (2 3 1) (3 1 %)■ Individual terms in this sum may contain infrared poles (or, after 
performing some of the integrations in complex plane, imaginary parts), but the expression as a whole is finite and 
real for e 7^ 0. 

We return now to the case of physical interest (to 2 = to 2 . = to 2 ; m 2 = 0) , and ignore all temperature-independent 
terms. We note, furthermore, that the contribution originating from the last term in Eq. (|A19(1 for (i, j, k) — (3, 1, 2), 
contains A[— (P3 + Pi) 2 ; to 2 ] + A[— (P3 — Pi) 2 ; to 2 ] which vanishes, given that P 2 + Pf = to 2 . The same is true for 
(i, j, k) — (2, 3, 1). This leaves us with 

Pf (to, /i) = [temperature-independent terms] + 

+ I h 11(0; to 2 , to 2 ) + 2P(to; fi) n(-m 2 ; to 2 , 0) + 

d d P i fd d p 2 n_(E 1 )n + (E 2 ) + n + (E 1 )n^(E 2 ) 1 



(2ir) d J (2ir) d 



SE\E 2 



Pl • p 2 - TO 2 - PiP 2 



d d Pl fd d p 2 n_(Pi)n_(P 2 ) +n+(P 1 )n+(P 2 ) 



1 



(2ir) d J (27r) d 



$EiE 2 



Pi • p 2 - to 2 + E 1 E 2 ' 



(A20) 
Pi • P2/IP1IIP2I, 



where we substituted pi — * — pi in the last term. We can still perform the integration over z 
leaving a rapidly convergent integral over |pi|, | p>2 1 ■ 

The final step is the expansion in e. The only temperature-dependent ultraviolet divergences are in the factorised 
terms on the second row in Eq. (|A20() . Adding together with the contributions from the other master integrals, as 
specified in Eqs. I|A5(I - (|A7|I . a straightforward computation reproduces the fermionic parts of Eqs. (| 1 1 11 — (f 1 3|) . 



APPENDIX B: FUNCTIONS DETERMINING THE MASS DEPENDENCE 

The functions that appear in Eqs. Ijll|l - (|13|l are defined as 



*i(v>A) 
Fa(v, A) 



1 



24tt 2 
1 



d.r 



dx 



X 


2 - 


x + y_ 





hv{\/x + y — A) + hp^x + y + A) 



X 


2 - 


x + y_ 





dx 



x + y 



ufWx + y - fi) + h F {^/x + y + A) 
UfWx + y - A) + n F (y/x + y + A) 



(i So ^ y/x~t + V\/x 2 + y 



(4tt) 4 



nsWx{ + y - p.)n F (y/x 2 + y + A) + n F (\fxi + y + fi)h F (y/xV+y - A) 



x In 



V x i + uV x 2 + y + y- \/xix 2 



V x i + V^/x 2 + y + y + yfxix 2 



+ 



np(\/ii + y - A)"f(v / ^2 + y - A) + %(V^7Ty + fj-)n F (y/x 2 + y + A) 

x In 



V x i + yV x z + y — y + \f x ix 2 ~ 



+ y\J x 2 + y -y- s/xTxi 



(Bl) 
(B2) 
(B3) 



(B4) 
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where 



rip (a;) 



1 

e x + 1 



(B5) 



These functions are related to the functions Jf, If and Hf defined in Eqs. I|A2(1 - I|A4() : the medium-modified part of Jf 
reads T A F\ for e = 0; the medium-modified part of If reads —T 2 F2 for e = 0; the medium-modified part of dlf/dr 
reads — Fs/^Air) 2 for e = 0; and the "non-factorizable" part of Hf (the last two terms in Eq. 
e = 0. The functions Pi, Fa, P3 are related by 



F 2 (y,A) = -2 



dy 



P 3 ( 2 /,A) = (47r) 



2 dF 2 (y,fi) 



c)y 



The functions introduced possess some solvable limiting values. For y 

Fi(0,A) 
F 2 (0,A) 
F 3 (0,A) 



7tt 2 




720 " 


24 


1 


A 2 


24 + 


8tt 2 


7T^ 


+ 2 7 



48tt 2 



23 



= In ■ 



16tt 2 



V V2 2tt/ ^V2 2vr 



201) rcad s T Ft for 
(B6) 

(B7) 
(B8) 
(B9) 



where 
setting y - 



denotes that the logarithmic divergence displayed on the right-hand side needs to be subtracted before 
* 0, and the function 23, which has the property 23(0) = 0, corresponds to the notation in Ref. [iof . 
The analytic expression in terms of ip{z) = T'(z)/T(z) comes from Ref. |2jj. For jj, ^ 0, the function F4 diverges 
logarithmically at small y, but as it is always multiplied by y, this behaviour has little interest in the present context. 
Inserting the values in Eqs. dHZJ-dHU into our expressions for aff, aff, aff, Eqs. (3.14), (3.15), (3.20) of Ref. Hj 
are reproduced. 

For y — > 00 and jl fixed, the functions Fi,F 2 ,F^ vanish as exp(— y/y), the function F4 as exp(— 2^/y), modulo a 
powerlikc prcfactor. An interesting limit is obtained, however, by setting y, /t simultaneously to infinity but keeping 
the ratio z = y/fi, 2 = m 2 /^ 2 fixed. This corresponds to setting the temperature to zero but keeping m, fj, finite. Then 



li m T 4 Ft ( —~ 

T^O V T 2 



limT 2 F 2 f^,a = 



0(1 
0(1 



Z) 8^ /2 ' 



limF,^ a = en- 



T-+0 



V T 2 

limT 2 F 4 f fflJ 

T->0 



V T 2 



^)-(/2 - W) 

z 



) = 0(1 -z) 



64tt 4 . 



-{w 4 



/ 2 2 ) , 



where 



f2 = VT 



z In 



(B10) 
(Bll) 
(B12) 
(B13) 

(B14) 



Inserting into our expressions for Qgj, Ojff , Eqs. (1) and (4) of Ref. [2g] are reproduced. 

Unfortunately the limit (1 — > 0, of most interest to us in this paper, does not render any of the functions analytically 
solvable, as far as we know. We show the results of numerical evaluations in Fig. |SJ 
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